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Scope of NASA-Ames Research Grant No. NAG 2-727 


Work under this agreement started on August 1, 1991. This final report is presented in 
two parts. The first part refers to computing the helicopter trim settings of periodic initial 
states and control inputs sequentially and in parallel. The second part refers to an 
exploratory study of a subspace iteration method as an alternative to the QR method for 
Floquet eigenanalysis. The papers resulting from this study are: 

1) Achar, N. S., and Gaonkar, G. H., "Helicopter Trim Analysis by Shooting and Finite 
Element Methods with Optimally Damped Newton Iterations", (In press) American 
Institute of Aeronautics and Astronautics Journal, 1993. 


2) Achar, N. S., and Gaonkar, G. H., "An Exploratory Study of a Subspace Iteration 
Method as an Alternative to the QR Method for Floquet Eigenanalysis", submitted for 
possible publication in the Journal of American Helicopter Society. 


PART - 1 


TRIM ANALYSIS 


G. H. Gaonkar 
N. S. Achar 



Helicopter Trim Analysis by Shooting and Finite Element 
Methods with Optimally Damped Newton Iterations 

Abstract 

Helicopter trim settings of periodic initial state and control inputs are investigated for 
convergence of Newton iteration in computing the settings sequentially and in parallel. The 
trim analysis uses a shooting method and a weak version of two temporal finite element 
methods with displacement formulation and with mixed formulation of displacements and 
momenta. These three methods broadly represent two main approaches of trim analysis: 
adaptation of initial-value and finite element boundary-value codes to periodic boundary 
conditions, particularly for unstable and marginally stable systems. In each method, both 
the sequential and in-parallel schemes are used and the resulting nonlinear algebraic 
equations are solved by damped Newton iteration with an optimally selected damping 
parameter. The impact of damped Newton iteration, including earlier-observed divergence 
problems in trim analysis, is demonstrated by the maximum condition number of the 
Jacobian matrices of the iterative scheme and by virtual elimination of divergence. The 
advantages of the in-parallel scheme over the conventional sequential scheme are also 

demonstrated. 


Notation 


a =Lift curve slope 

c =Control-input vector 

C d = Resultant profile drag force in the plane of the rotor disk opposite 

to the flight direction 
C d0 =Profile drag coefficient 
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=Rolling moment coefficient 

=Pitching moment coefficient 

=Thrust coefficient 

=Blade pitching-moment coefficient 

=Weight coefficient of the helicopter 

=Equivalent flat plate area of parasite drag 

=Aerodynamic moment per unit length of the blade in flap and lag 
directions, respectively 
=Objective function to be minimized 
=Hamiltonian 
=Identity matrix 

=Jacobian or nondimensional torsional inertia 
=Lagrangian 

=Generalized momentum (L^) 

=Flap natural frequency (rotating) 

= Generalized coordinate 
=Nonconservative force 
=Rotor radius 

=State vector y augmented with control-input vector c 
=Initial and final times 
=Kinetic energy 

= Flight speed or Potential energy 
=Work done by Q 
=State vector 

=y(t) with initial condition y(0) = yo 
= Shaft tilt angle 
=Flap response 
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=Lag response 
= Azimuth angle 

=Collective, longitudinal cyclic and lateral cyclic pitch angles, 
respectively 

=Pitch angle ( = 0 o + 0 C cos^ + 9s sin ip ) 

=Advance ratio ( = V cosa s /^R ) 

=Dimensionless flight speed ( = V/flR ) 

=Dimensionless nonrotating flap, lag and torsional natural frequencies 
=Rotor speed 

=Lock number (blade inertia parameter) 

=Newton damping parameter 
=Inflow 

=Newton direction (see Eqs. 4 and 5) 

=Rotor solidity 
=Transpose of [] 

=Vector norm 
=Gradient of g 
=Time derivative of ( ) 

=Partial derivative of ( ) w.r.t. q, similarly subscripts p, q, q, c and y 
indicate partial differentiation. 

Introduction 


The helicopter trim settings comprise control inputs for required flight conditions and 
the corresponding initial conditions for periodic response. They are prerequisite for stability 
and vibration studies. The control inputs are specified indirectly so as to satisfy flight 
conditions of prescribed thrust levels, rolling and pitching moments etc. In addition to the 
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nonlinearity of the system and control-input equations, the control inputs appear not only 
in the system damping and stiffness matrices but in the input matrix as well, and must be 
found concomitantly with the periodic response. The prediction of trim settings has been 
vigorously pursued since the 1980s and still is a demanding exercise because of divergence 
of iterative schemes and excessive machine time (Refs. 1—6). 

Particularly, for marginally stable and unstable systems, the shooting method is 
increasingly used (e.g., 2GCHAS, Ref. 7); however, much recent research has been centered 
on temporal finite element methods of different versions, such as displacement, mixed and 
bilinear formulations, with further classifications in each of these formulations involving h, 
p and hp versions (Refs. 4, 5, 8, 9). No matter which method is used, computation of trim 
settings leads to nonlinear algebraic and transcendental equations, whose solution at 
present cannot be based on solid theory (Ref. 10). In fact, the computational difficulties of 
these equations virtually preclude the translation of several trim analysis methods into 
robust algorithms with global or reasonably qualified convergence characteristics. Little 
information is available on the nature of such difficulties or on ways to quantify and 
alleviate them. Newton’s method is the most widely used and perhaps the best method of 
solving nonlinear equations (Ref. 10). But while it guarantees quadratic convergence (the 
number of significant or accurate digits doubles after each iteration), it guarantees only 
local convergence and is sensitive to the initial guesses or starting values. And even with 
good starting values, the method can exhibit erratic divergence due to numerical corruption 
(Ref. 1). 

The present investigation covers this divergence problem with respect to the shooting 
and two temporal finite element methods, which typify broadly two classes of methods: 
adaptation of time marching methods of initial-value problems and finite element methods 
of boundary-value problems to periodic boun da ry conditions. It also covers the in-parallel 
scheme or the simultaneous computation of initial conditions and control inputs vis-a-vis 
the sequential scheme (Ref. 2), in which the iterations for the initial conditions and control 
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inputs are carried out as two separate computational blocks, one following the other (Refs. 
1,2). It must be emphasized that the divergence problem is not peculiar to the in-parallel 
scheme. In fact it is as much a part of the sequential scheme. Moreover, the bulk of the 
earlier trim analysis investigation uses sequential computation. An exception is Ref. 1, 
which found appreciable machine time saving through the in-parallel scheme in the 
shooting method. However, that finding is masked by erratic divergence of Newton 

iteration. 

Given this background, the present investigation is noteworthy in the following respects. 

1. The Newton damping parameter is examined concerning both its selection (with a 
rational basis of minimizing an objective function) and its role in alleviating the sensitivity 
of Newton iteration to the starting values in the solution of trim settings of initial state 

and control inputs. 

2. The computational reliability of the Newton iteration without and with optimal 
damping is quantified by the condition number of the Jacobian matrix, which also explains 

rationally the earlier-observed divergence problems (Ref. 1). 

3. Concerning divergence and machine time, a comprehensive comparison of the 
sequential and in-parallel schemes is provided; each scheme is treated with Newton 
iteration both without and with damping. This exercise includes three trim analysis 
methods, representing two main approaches of trim analysis, particularly in stability 
investigations. 


Damped Newton method 


The method retains the highly attractive features of the original Newton method (e.g., 
quadratic convergence) and yet almost global convergence (Ref. 10). We consider the 
solution of n nonlinear equations 

f i(s i, s 2 j •••i s n) = 0; i = l,2,...,n 


( 1 ) 
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or, in equivalent form, 

f(s) = 0 ( 2 ) 

for which the Jacobian matrix is given by 

J ( S ) = 3T = ; ^ -» n ^ 

The algorithm begins with the "improved" solution 

s m+i = s m + (4) 

where m is the iteration counter, A is the "optimal" damping factor and £ is the solution of 


the linear system, 

f = — J(s)" 1 f(s) (5) 

The terms improved, and optimal are qualified in the absence of a solution s* such that f(s*) 
= 0 and of optimality conditions to determine A. The theory of unconstrained minimization 
and weak line search (Refs. 10,11) provides a rational basis of quantifying these two terms 
and solving for the damping factor. We bypass the mathematical details and include 
instead a brief account of the method, following Ascher et al. (Ref. 10). 

gm+i is an improvement over s m in the sense of minimizing an associated objective 
function g(s m + A£) monotonically, where 

g(s) = 0.5j i f i (s) 2 (6) 

The objective function has the property that g(s) > 0 and g(s*) = 0 when f(s*) = 0. Thus, 
the minimum of g(s) provides the solution. Moreover, the Newton direction is a descent 

direction; that is, for the gradient Vg, we have 

£tv g = _[j-if]t[jtf] (7a) 

= - f(s) 2 = -2g < 0 (7b) 

Expanding g(s+A£), we obtain 

g( s m + Af) = g(s“) + Aevg(sm) + 0{ A2 | ^ | 2 ) < g(s“>) (8) 

where Vg, which is equal to J l f, shows that minimization is sought in the Newton’s 
direction. Thus s m+1 is an improvement over s m in the sense that 

g( s m*i) = g( s m + A£) < g(s m ) 


( 9 ) 
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which says that the solution s* is keyed to the generation of monotonically decreasing 


values of g(s). 

We define from Eq. (8) 


g(s m + Af) - g(s“>) 

<p(\) 

A^Vg(s m ) 


The algorithm (Refs. 10,11) makes use of a fixed parameter a such that 


0 < a < 0.5 
a < <^(A) < (1— o') 


( 10 ) 


(Ha) 

(lib) 


From Eqs. (8) and (11), we have 

(1—2 A (l-a)) g(s“) < g(s m+1 ) < (1— 2Aa) g(s m ) (12) 

The role of o is sketched in Fig. 1, which shows the inherent predictor-corrector structure 
of the algorithm (Ref. 10). For the mth iteration counter, we approximate the scalar 

objective function by a quadratic: 

g(s m + A£) a aA 2 + bA + c = ^(A) (13) 

where the three unknown constants — a,b and c are determined from the following 


conditions: 

ip(0) = g(s m ) (14a-) 

V*Am) = g(s m + Am£) ( 14b ) 

^(°) = g I ( sm + AOI o =^Vg(s“) (14c) 

The criterion that ^ & minimum with respect to A can be expressed as 

-AJ ^'(0) Am (15) 

A = — <- 

2('0(Am) — ^(0) — A m ^ / (0)) 2(1— cr) 

The preceding algorithm has been found to work generally well (Ref. 10); direct 
computational experience in computing trim settings supports this as well. However, there 
are cases of the algorithm breaking down when the objective function fails to decrease 
monotonically. That is, when g(s) has local minima and/or singularities, Eq. (13) may not 
be a good approximation to g(s). Hence, the algorithm fails to compute A satisfying Eq. 
(9). This problem has been alleviated as follows. At the end of every iteration, before 
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computing the damping parameter A, s^i = s m + Af is computed with A = 1.0. Then, if 
any of the control inputs in s m+1 exceeds physically realistic values (say shaft tilt <2s,max — 
20°), an upper limit for A, A., is chosen such that the control inputs in s* + i are physically 
realistic. Then, the algorithm proceeds to find optimum A in the range 0 to A m as before. 
Again, in extreme cases if the algorithm fails to find A satisfying monotonicity Eq. (9), A m 
itself is chosen as the optimum A. Though the monotonic decrease in g(s) is not guaranteed 
with A = A m , the algorithm converges to the trim settings, which are elaborated on later 
with the help of numerical results. 

Condition number of J 


The relative error in the solution of trim equations by Newton iteration can be bound by 
utilization of the condition number of the Jacobian matrix, cond(J), which also quantifies 
the robustness of the Newton direction; see Eq. (5). To provide an improved appreciation 


of the role of cond(J), we emphasize that the actual computation of trim equations does not 
follow Eq. (5); it follows a numerically perturbed equation: 

[J + <5J] {£ +8£} = {f + <5f} ( 16a ) 

The following inequality (Ref. 12) 


m 

lien 


■ < cond(J) 


\m\ 


Pll 


■ + 


11**11 


(16b) 


shows that cond(J) represents the maximum possible magnification of the sum of relative 
errors in J and f. Thus, the higher the value of cond(J), the greater the sensitivity of Eq. 


(16a) to computational perturbations and, consequently, the less well-conditioned is the 
computational problem of finding the control inputs and periodic initial state. 


With the definition 


6 = max { l/cond(J) } 


and with Eqs. (3) - (7), it can be shown that (Ref. 10) 

6 |Vg| U| = <5| J 1 f| |J-‘f| < |f| 2 = -Vg^ 


(17) 


( 18 ) 
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Equations (5) and (18) show that, with increasing value of cond (J), the overall 
conditioning of the Newton direction decreases. That is, the product V g‘f decreases 
although the correction |f[ is not small and g(s) does not decrease rapidly along the 
Newton direction f. Schematically stated, lines ab and ac tend to merge with line ad in 
Fig. 1. 


Trim formulation 


We include a brief account of the shooting method (Ref. 1) and the weak version of a 
temporal finite element method with mixed formulation (Ref. 8) of displacements and 
momenta. The algorithmic details of the temporal finite element method with displacement 
formulation run similar to those of mixed formulations and are omitted here; for details see 
Ref. 4. This facilitates appreciation of the algorithmic aspects of sequential and in-parallel 
schemes of Newton iteration in the trim analysis by the shooting and finite element 
methods. For convenience, the latter two finite element methods of mixed and 
displacement versions are, respectively, represented as FEM— M and FEM— D. In the three 
methods, the algorithm follows the in-parallel scheme or the simultaneous computation of 
initial conditions and control inputs. The straightforward adaptation to the sequential 
scheme is not spelled out explicitly. 

Shooting method 

In trim analysis, equations of motion in state variable form 

y = G(y(t),c) (21) 

satisfy the unknown periodic initial state yo, that is, 

y(27r;y 0 ) - yo = 0 (22) 

Further, the unknown control inputs c should be determined such that the desired flight 


conditions 


f(y,c) = o 


(23) 



10 


are satisfied. Equations (22) and (23) comprise the nonlinear algebraic trim equations, 
which are symbolically represented as Eq. (2), where s = |_y,cj 1 is the augmented vector of 
trim settings. These equations are solved using Newton iteration to get the trim settings. 

FEM-M 

The Hamilton’s Law of Varying Action or HLVA of a system can be represented as 

tf t f 

6 f( L + W) dt - (L • 6q)\ =0 (24) 


to 


to 


(25) 


in which the system is represented in configuration space in terms of the generalized 
coordinates. In the mixed formulation, we first represent Eq. (24) in terms of the 
coordinates from the phase space using the Hamiltonian of the system, such that both 
generalized displacement q and generalized momentum p become primary variables. By 
varying the Hamiltonian of the system, 

H = H (q,p,t) = p.q-L (q, q, t) 

we obtain 

£H = (5p.q +p.<5q — <5L(q,q,t) 

Hence, 

<5L = <5p.q + p.<5q - <5H(q,p,t) 

= tfp.q +p.<5q - [ H q .<5q + H p .tfp ] 

Substituting for 6L in Eq. (24), we have 


tf 


f ( £pq +P-£q - [ H q .(5q + H p .5p - Q <5q ] ) dt = (p.£q) | 


(26) 


to 


In the above equation, q occurs only in the first term. Hence, integrating the first term by 

parts, q can be eliminated from the above equation to get 

tf tf 

/ ( ~ <5p q +P <5q - [ Hq.tfq + H p .£p - Q <5q ] ) dt = [ (p.<5q) - (q.£p)] | (27) 


to 


to 


With the definition 
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and 
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'6q 
8 P 


H = 



l J 
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P 

-q. 


( 28 ) 


Eq. (27) can be expressed in vectorial form: 

f ( 6y.Ny - 6y.il ) dt = [6y.B] | (29) 

t to 

to 

When H, defined in Eq. (28), is nonlinear in y, it can be linearized about a steady state 
value y as 

H = H(y) + H(y) y Ay (30) 

where 


H(y) = 


H q -Q 

H n 


, H(y) y = 


Hqq - Qc 


H qp - Qp 


H 


pq 


H 


pp 


(31) 


Substituting Eq. (30) in Eq. (29), we obtain 


tf 

f [ <fy-N (y + Ay ) - Sy H(y) - 6y [H(y) y ] Ay ] dt = (ffy.B)| (32) 

to ° 


Next, the time interval [t 0 ,tf] is discretized into m smaller segments; i.e., t 0 = ti < t 2 < --ti 
• •< t m+ i = tf. In each of these temporal elements, the generalized coordinate q e and the 
momentum p e can be expressed in terms of some appropriate shape functions as follows. 

Since the derivatives of q and p do not occur in Eq. (27), constant values with discrete 
end-momenta and displacements will satisfy the completeness requirements (Ref. 8); i.e., 


we can have 
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ye = 


'q 

p 


Qi 

Pi 


if ti < t < t i+1 


fqil 


ye = 


Pi 


if t = ti ; y e = 


'quf 

Pi + 1 


if t = ti + j 


(33) 


However, the virtual displacement and the virtual momentum 
differentiable functions (Ref. 8). Hence, we choose 


<fye = 


'«q 

r 

6p 

- 


(1-r) 0 t 0 

0 (1-r) 0 t 


' £qi 

1 . foi 

•* ^qi+i 
^ ^Pi + l 


= M 6y e 


require piecewise 


(34) 


with 

(t - ti) 

r = (ti + r^ti) 


Substituting for Ay and 5y in Eq. (32) from Eqs. (33) and (34), respectively, for the ith 

time element, we have 
twi 

6y e [ f [Mt.N (y e + Ay e ) - M‘ H(y) - H(y) y Ay e ] dt - (Mt B) | ] = 0 (35) 

t i 

Defining 

t i+i ti+i 

Fe = / [ Mt H(y) - Mt.N y e ] dt, K e = f [M l H(y) y - Mt.N ] dt 

ti ti 


and 


t i + l 

Ge = (Mt B) I - 
t i 


-Pi 

qi . 

Pi + 1 
— qi+i 


Eq. (35) can be expressed as 

«5y e [ F e + K e Ay e + G e ] = 0 (36) 

Next, all the m elemental Eqs. (36) are generated and assembled to get the global, 


linearized variational statement 
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6y [ F+ K Ay+ G] — 0 

where K, F, G, y, <5y are the global stiffness matrix, force vector, momentum vector, 
displacement vector and the virtual displacement vector, respectively. The elements of the 

vectors G, y and 8y are given below. 

G = [ — pi, qi, 0, ... 0, 0, p m+ i, — qm + i J 1 (38a) 

y = L Ql> Pl> Qa> P 2 ,-", Qm, Pm J 1 ( 38b ) 

6y = [ <5qi, tfpi, <5q2, 8p 2 &lnui» £pm + i J 1 ( 38c ) 

Then, applying the periodic boundary conditions to Eq. (37) (i.e., qi = q m+ i and pi = Pm+i), 

we get 

F + K Ay = 0 ( 39 ) 

where F and K are F and K, respectively, rearranged after applying the boundary 
conditions. Next, Eq. (39) is solved using Newton iteration until the series y = S Ayj 
converges, both F and K being updated at the end of each iteration. 

Parallel trim method 

For the parallel solution of trim settings, we reformulate the F and K of Eq. (39) as 
follows. Since H, defined in Eq. (28), is nonlinear in both y and c, it can be linearized about 
some mean position y and C. (Note: The Hamiltonian H is nonlinear in q and p, and the 

generalized force Q is nonlinear both in p, q and c.) That is, 

H = [ H(y,C) + H(y,C) y Ay + H(y,C) c Ac ] (40) 

where 

Hqc — Qc 
H(yA- = H pc 

and H(y,C) y = H(y) y ; see Eq. (31). Substituting Eq. (40) in Eq. (29), we get the 

variational statement for each element as 
t i + l 

6y e [ f [MbN (y e + Ay e ) - Mt H(y) - H(y) y Ay e - MtH(y,C) c Ac] dt 

t- 

M + l 

- (Mt B) | ] = 0 

t i 


( 41 ) 
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Defining 

1 i + l 

F e = / [ Mt H(y,C y ) -M‘.N y e ] dt, 
t i 


t i + l 

K e = f [M l H(y,C) y -M‘.N ] dt 
t i 


and 



-Pi 
Qi 
Pi + l 

^ — Qi4l 


(42) 


Eq. (41) can be represented as 

<5y e [ F e + K e Ay e + Kec Ac + G e ] = 0 (43) 

Next, all the elemental Eqs. (43) are generated and assembled to get the global, linearized 
variational statement 


<5y [ F+ K Ay + K c Ac + G] = 0 (44) 

Further, Eq. (23) can be linearized as 

f(y,c) = f(y>c) + f(y,e)y Ay + f(y,c) c Ac 
Now, Eqs. (44) and (45) are combined and the boundary condition is applied to get the 


= 0 


(45) 


augmented force vector and stiffness matrix, F and K of Eq. (39), respectively. Then, Eq. 
(39) is solved iteratively, until the augmented vector s = EAsi converges where Asi = [Ayi, 


AciJ t. 


Model description 

For computational purposes, flap— lag and flap— lag— torsion models are selected, both the 
models are based on quasisteady aerodynamics and rigid body mode representation. 
However, for the simplicity of illustrating the algorithmic and computational aspects, 
model description and much of the discussion of numerical results are for the flap-lag 
model, which was also treated in Ref. 1 by the shooting method. We begin with the state 
vector y(t) with four components comprising the flap angle /? and lag angle ( and their 
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rates 0 and y(t) = l /J(t), 0(t), <(t), C(t) J >. Part of the trim analysis is to compute the 
periodic initial state yo such that 

y(27r;y 0 ) - yo = 0 ( 46 ) 

The remaining part of the trim analysis is to compute the four parameters — three pitch 
angles 0 0 , 0 S and 0 C and shaft tilt a s — to satisfy the following four trim equations of force 


and moment balance: 


C t cos(a s ) + Cd sin (a s ) = C w 
Ct sin(a s ) - Cd cos(a 5 ) = 0.5 2 f 
Ci = 0 
C„, = 0 


(47) 

(48) 

(49) 

(50) 


The solution of four initial-condition equations, typified by Eq. (46) and four trim Eqs.(47) 
— (50), constitutes the trim analysis. The n ondim ensional thrust Ct and horizontal force Cd 
depend on the total blade root shear forces: F p , normal to the blade in the flap direction; 
F<., in the lag direction, and F r , in the outward radial direction. These shear forces are 

given by 


F p = -1.5 - 1.5 cos{0) sin(/3)(l-K) 2 + / F p dr 


(51) 


F^ = —1.5 £ cos 2 (/3) + 3 sin(/?)cos(/?) (1+C) ft + cos 0 f F C, (^) 

0 

F r = 1.5 [ (1+C) 2 cos 2 (/?) + P 2 ] (53) 

Then, C t and Cd, in Eqs. (47) - (48), and the rolling moment and the pitching moment 
coefficients, Ci and C ra , respectively, in Eqs. (49) - (50), are given by 


27T 


CiL 

era 


er a 

J. 

7 27t 


= — / [ p p cos(0) + F r cos(/?)]dV> 


(54) 


7 
27T 


f [Fp sin(/3)cos(^+0 + F<- sin(^+C) + F r cos(/?)cos(^+C)] d^ (55) 

n 

( 56 ) 


Ci _ 
era 


(P,?-l) . 

= f PMi>+OW 
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2tt 

J* P cos(^H-C) dV' 


The final component of the trim analysis concerns the uniform total inflow X\: 


A i = fi tan(a s ) + C t / [2 + M )] 


which is solved iteratively in combination with Ct and a s . 


(57) 


(58) 


Trim analysis results 


Trim analysis uses the shooting method and the two finite element methods, FEM M 
and FEM-D. The control inputs and initial conditions are computed simultaneously as one 
block in the in-parallel scheme; they are computed sequentially as two separate blocks in 
the sequential scheme. Unless stated otherwise, the rigid flap-lag model with in-parallel 
scheme is used; nonlinear equations are solved by conventional Newton iteration with no 
damping and the following baseline values are assumed: 7 = 5, o>p = 0.57, = 1.4, a - 

0.05, a = 6.28, C w = 0.01, C do = 0.01, f = 0.01, 0 < Ji < 0.7. The computations are 
performed on a VAX 6320. The sparse matrices obtained by the finite element methods are 
solved using the NAG subroutine (F01BRF), wheih considers the matrix sparsity. 

In the finite element methods, it is first necessary to arrive at the number of elements, 
NEL, needed for a priori specified level of tolerance in the solution of periodic response. 
This tolerance is further substantiated on the basis of a relative error norm criterion with 
the shooting-method results (Refs. 1,2) as reference values. For that purpose, a typical 
flight speed, /x = 0.4 is chosen. As shown in Fig. 2a for a typical initial state P(0) t all the 
periodic initial conditions, y = [ /?(0), /?(0), ((0), C(0)J converge as the number of elements 
increases. In particular, it was found that for asymptotic convergence of all the four 
components of y, we need at least NEL = 12 for FEM— D and NEL = 16 for FEM— M. 
Concerning the minor differences between results from the three trim methods, the relative 
error norm in y = [ P, P, C (J \ obtained by the respective FEM, is defined as 
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Relative error norm = 




5j 1 1 yshoo t yfem || j 
S|| yshootll i 


], i = 1,...,NEL 


As seen from Fig. 2b, as the number of elements increases, the relative error norm 
decreases rapidly; for NEL = 12 for FEM-D and for NEL = 16 for FEM-M, the relative 
error norm is less than 0.01. For the above NEL values selected on the basis of the results 
at Ji = 0.4, it was also verified that right up to fi = 0.7, the relative error norm is less than 
0.01. Thus, in summary, NEL values of 12 and 16 for FEM-D and FEM-M, respectively, 
guarantee 1) asymptotic convergence for n = 0.4 and 2) a relative error norm of less than 
0.01 for 0.0 < Jt < 0.7. Spot checks for other values of /i show that asymptotic convergence 
with these NEL values holds for 0.0 < Ji < 0.7 as well. Overall, with these NEL values, the 
control inputs 8 0 , 0 S , 0 C and a s , and the periodic responses /?, ft, £ and £, agree with the 
shooting method results. Hence, these NEL values are chosen in all the subsequent 
numerical results with the flap-lag model. ( This agreement is further elaborated later on 


for the flap-lag-torsion model.) 

In Figs. 3-5, the machine times taken by the sequential and in-parallel schemes in the 
three trim analysis methods are presented. Given the sensitivity of Newton iteration to 
starting values, the results for each of the methods are presented for two sets of starting 
values. In part (a), we use the "exact" solution of the preceding flight speed Ji as the 
starting value; the exact solution is taken as the one obtained by continuation approach 
with Ji as a continuation parameter and A ft - 0.05 as the continuation step size. For 
example, starting values, say at Ji = 0.3, are given by the solution at the preceding value of 
Ji = 0.25. For a given flight speed, the cumulative machine times taken starting from n = 0 
are shown in parts (a) of Figs. 3-5. Though prohibitively costly, this approach to the 
starting values provides a rational basis of providing the ‘best’ starting values. The other 
extreme is to begin with zero starting values, perhaps the most demanding starting values 
for the iteration. This is done in part (b) of Tigs. 3-5; owing to the divergence problem of 
Newton iteration with zero starting values in both the sequential and in-parallel schemes, 


18 


the results are limited to /z < 0.3. Figures 3—5 show that the in-parallel scheme is more 
economical than the sequential scheme. This saving is observed for both sets of starting 
values (continuation with A Ji = 0.05 and zero starting values), showing that the in-parallel 
scheme is preferable regardless of the starting values. 

The preceding results are based on Newton iteration with no damping, and mention is 
made of the divergence of the iteration with zero starting values. The impact of damped 
Newton iteration on divergence and related issues are pursued in Figs. 6-9; an important 
observation is that the damped Newton method did not encounter divergence. Figures 6 
and 7 show the mechanism of divergence relative to iteration counter and the starting 
values. This is followed by Figs. 8 and 9, which show a means of quantifying and 
understanding divergence as well as computing the Newton damping parameter, which 
virtually eliminates divergence. 

Figure 6 shows the iteration counter versus flight speed. While iteration without 
damping experiences divergence for approximately /z = 0.45, damped iteration does not 
diverge although its iteration counter increases with increasing /z. That the Newton 
damping parameter makes the iteration more controlled and "gives less room for erratic 
behavior" (Ref. 10) is well borne out in Fig. 6. Also, for /z > 0.6 or so, the iteration 
counter in the damped iteration rapidly grows in the shooting method, indicating poor 
convergence. By comparison the damped iteration in FEM— M and FEM— D is remarkably 
smooth; the iteration counter and its growth with increasing Ji are much less rapid and it 
hardly exceeds 15 for the complete sweep of 0.0 < n < 0.7. As seen from Fig. 6, the 

FEM-M and FEM— D have fast convergence up to (i = 0.7 while the shooting method 
exhibits slow convergence for Ji > 0.65. Summarizing, we observe that the FEM-M and 
FEM— D are much better regarding iteration counter or speed of convergence; at /z = 0.7 for 
example, iteration counter for FEM— M and FEM— D is about 15 whereas it is about 55 for 
the shooting method, nearly four times higher. 

The divergence problem of the Newton iteration due to the starting values is further 
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pursued in Fig. 7 in th ep-p plane, where p is related to the starting values as 

starting value = p x exact solution 

Here, p = 0.0 and p = 1.0 imply that all the starting values are zero and exact solution, 
respectively. (The exact solution at any flight speed is that solution obtained by 
continuation approach with Ap = 0.05). This is done as a means of quantifying the 
sensitivity to possible extremes of starting values and of connecting these results with an 
earlier investigation (Ref. 1). The other possibility of large initial values (p > 1) is not 
exercised separately. The divergence boundary of Fig. 7 corresponding to the shooting 
method is similar to that of Ref. 1; the minor differences are due to the sensitivity of these 
boundaries to discretization in p and p. Even with somewhat improved starting values, say 
p — 0.1 compared with p = 0.0, erratic behavior of the boundaries merits mention. The 
damped Newton iteration reduces the erratic behavior in general; indeed in Fig. 7 it 
converges everywhere. With poor starting values (low values of p, say less than 0.5), 
divergence with Newton iteration is not unexpected since it guarantees only local 
convergence and is sensitive to starting values anyway. What is unexpected is divergence in 
the shooting method for p as high as 0.8. Concerning divergence boundary, shooting 
method is affected more than FEM-M and FEM-D. This supports the iteration counter 
results in Fig. 6. 

The unexpected divergence with Newton iteration and the absence of divergence with 
damped Newton iteration is further investigated in Fig. 8 on the basis of the maximum 
condition number of the Jacobian matrices of the iterations; see Eqs. (16a) and (16b), 
respectively. The results are for p = 0.0, zero starting values. As seen from Fig. 8, the onset 
of divergence in Fig. 6 for p = 0.50 in the shooting method and p = 0.45 in FEM-M and 
FEM— D is accompanied by rapid increase in the Jacobian matrix condition numbers for 
corresponding values of p. By comparison these condition numbers essentially remain 
constant in all the three trim analysis methods with damped Newton iteration. This shows 
that damped Newton iteration significantly improves the overall conditioning of the 
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iteration. 

Figure 9a shows, for Ji = 0.3, the monotonic decrease of the objective function, see Eq. 
(6). We emphasize that g(s) = 0 gives the damping parameter A. It is seen that the 
objective function decay is monotonic on expected lines in all the three trim analysis 
methods and that the minimization is fairly rapid. That is, the iteration counter hardly 
exceeding 7 for g(s) = 10' 11 for each of the three trim analysis methods. Figure 9b shows 
the initial nonmonotonic decrease of the objective function. In such cases the algorithm 
follows the remedial measures referred to earlier in conjunction with Eq. (9). Two 
distinguishing features of Fig. 9b merit mention. First, even in cases when A cannot be 
computed on the basis of monotonic decrease of the objective function, the algorithm still 
converges to the trim settings with as few as 7 iterations in all the three methods. Second, 
despite a couple of initial iterations involving nonmonotonic decrease, it later, converges 
with a rapid monotonic decrease. 

Thus far, the influence of the damped Newton iteration on the divergence problem of the 
in-parallel scheme is addressed. We now conclude with a brief discussion of its influence on 
the sequential scheme. In this scheme, unlike the in-parallel scheme, there are two sets of 
starting values; one for the response loop and one for the control loop. Without damping, 
the divergence boundary for the shooting method with the sequential scheme is given in 
Fig. 10; for the other two methods they remain nearly the same and are not shown. For all 
the three methods, the starting values for the response loop are assumed zero and those for 
the control loop are chosen using the starting value parameter p : 

Starting value = px exact control settings 

where the exact control settings are those obtained by the continuation approach with A p, 
= 0.05. In this scheme, it is observed that the divergence is mainly due to the sensitivity of 
the response loop to the physically unrealistic values of the control inputs, estimated in the 
control loop. When the control-input computations are prevented from generating 
unrealistically large values either by damping or by a priori stipulation, the response is 
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always found to converge. That is, damping the control loop is more effective than 
damping the response loop, and, therefore, in this illustrative example, the damped Newton 
iteration is implemented only in the control loop. After introducing the damped Newton 
iteration, in both versions of FEM, the sequential scheme converges in the entire p - Ji 
plane right up to Ji < 0.7. However, in the shooting method, i.e. Fig. 10, though there is an 
increase in the region of convergence because of damping the Newton iteration, the 
iteration was found to cycle for p beyond 0.55; i.e., after some iterations, the algorithm 
computes a solution that it has already computed in one of the previous iterations and 
hence the algorithm enters into an infinite loop without converging. This cycling 
phenomenon is found to be independent of the ill-conditioning or near-singularity of the 
Jacobian as quantified by the condition number of the Jacobian. This problem is involved 
and not well-understood, (Ref. 10) and there seems to be no known method to treat it 
effectively. 

The preceding investigation based on the flap-lag model results is further verified on the 
basis of flap-lag-torsion model results. The rotor parameters are identical to those of the 
flap-lag model; the additional torsional parameters are = 3.0, J = 0.002 and = 
-0.02. Overall, the results are nearly identical to those of the flap-lag model results, 
specifically these results refer to damped Newton iteration (e.g., divergence boundary, 
condition number and monotonic decrease of objective function) and to sequential scheme 
vis-a-vis in-parallel scheme. In Fig. 11, we select for illustration the torsional response rate 
y>(t) among the six components of the state vector y = [/?, ft, (, (, tp, <^J 1 an d the shaft-tilt 
angle a s among the four components of the control-input vector c = |_0o> @s> ^sj • The 

results from the three methods agree. Finally, as in Figs. 3-5, Fig. 12 shows the machine 
time saving with the in-parallel scheme in the FEM — D method; similar trends (not shown) 
are exhibited by the FEM-M and shooting methods with the in-parallel scheme. The only 
major difference between the flap-lag and ' flap-lag-torsion results concern NEL or the 
number of elements needed for asymptotic convergence of the periodic response as shown in 
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Fig. 2. For the latter case, because of highly oscillatory nature of the torsional response, 
NEL values as high as 60 and 90 are needed for FEM-D and FEM-M, respectively (Ref. 
14). 


Concluding remarks 


The three trim analysis methods are investigated with an optimally selected damping 
parameter A; A = 1 refers to conventional Newton iteration. With each trim analysis 
method, both the sequential and in-parallel schemes of computations of initial conditions 
and control inputs are exercised. The computational efficiency is described on the basis of 
both machine time and convergence characteristics, which are quantified by the maximum 
condition number of the Jacobian matrices of the Newton iteration. The iteration counter 
and its growth with increasing flight speed and divergence boundaries correlate with this 
quantification. That investigation demonstrates the feasibility of using an optimally 
selected Newton damping parameter in the in-parallel scheme to improve the 
computational efficiency of the trim analysis. It also shows the following. 

1) In the three trim analysis methods with both the sequential and in-parallel schemes, 
the optimally selected damping parameter virtually eliminates divergence up to flight speed 
Jl = 0.7 except for a small region beyond fj. > 0.55 in shooting method with the sequential 

scheme. 

2) The in-parallel scheme takes much less machine time compared to the sequential 
scheme with comparable convergence characteristics. 

3) At very high advance ratios (for Ji > 0.6 or so), the shooting method shows slow 
convergence in that the iteration counter and the machine time increase rapidly. By 
comparison, FEM-M and FEM-D show fast convergence for the entire range of 0 < n < 
0.7. 

4) The cycling phenomenon observed at /x > 0.55 in the shooting method with the 
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sequential scheme and with damped Newton iteration merits further research. 

The preceding investigation is restricted to a generic Newton-Raphson iteration and does 
not consider a wide class of related methods such as quasi Newton and other globally 
convergent methods. Nevertheless it should provide a reference point for using and 
comparatively assessing such methods in helicopter trim analysis. 
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Fig- 1 An acceptable range of the damping 
parameter A for the damped Newton iteration 
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Fig. 2 Convergence of Periodic initial-condition 
in FEM-D and FEM-M with increasing number of 
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Fig. 3 Machine time comparison in the shooting 
method with sequential and in-parallel schemes 
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Fig. 4 Machine time comparison in the FEM-D 
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Fig. 5 Machine time comparison in the FEM-M 
with sequential and in-parallel schemes 
( sequential, in-parallel) 








Fig. 6 Iteration counter comparison for the 
three trim analysis methods with zero sta- 
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Fig. 11 Comparison of the torsional response 
and the shaft tilt angle from the three trim 

analysis methods ( Shooting, FEM 
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FLOQUET EIGENANALYSIS 


G. H. Gaonkar 
N. S. Achar 


An Exploratory Study of a Subspace Iteration Method as an 
Alternative to the QR Method for Floquet Eigenanalysis 

Abstract 

Floquet eigenanalysis requires a few dominant eigenvalues of the Floquet transition 
matrix (FTM). Although the QR method is used almost exclusively, it is expensive for 
such partial eigenanalysis; the operation counts and, thereby, the approximate machine 
time grow cubically with the matrix order. Accordingly, for Floquet eigenanalysis, the 
Arnoldi-Saad method, a subspace iteration method, is investigated as an alternative to the 
QR method. The two methods are compared for machine-time efficiency and computational 
reliability, which is quantified by the condition numbers of the required eigenvalues and 
the residual errors of the corresponding eigenpairs. The Arnoldi-Saad method takes much 
less machine time than the QR method with comparable computational reliability and 
offers promise for large-scale Floquet eigenanalysis (say, FTM order > 100). 

Introduction 

Floquet theory is the primary mathematical tool in the investigation of rotorcraft 
stability. Its application involves computation of the following: 1) the trim settings of 
initial state and control inputs for periodic response satisfying the flight conditions; 2) the 
FTM about that trim response, and 3) a few largest eigenvalues of the unsymmetric FTM, 
which often is a byproduct of the trim analysis. This study examines the computational 
aspects of the third item. 

The information base on generating the FTM and its eigenanalysis is limited to relatively 
small-order systems (say, system or FTM order smaller than 100) and the generic QR 
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method is used almost exclusively for eigenanalysis (Ref. 1). Though the machine time 
needed to generate the FTM is sensitive to several factors, such as system nonlinearity, and 
defies generalization, it far exceeds the machine time for the eigenanalysis of the FTM. 
Therefore, for small-order systems, the search for a viable alternative to the QR method is 
not an issue. But for large-order systems, this search merits further investigation for two 
related reasons. First, the QR method is the recommended method for a complete 
eigenanalysis while Floquet eigenanalysis requires only a few most dominant eigenvalues or 
partial eigenanalysis. Second, for large systems, the QR method is expensive since the 
operation counts grow cubically with the matrix order (Refs. 2,3). Therefore, the machine 
time also grows approximately cubically with the matrix order. Further, the recent 
developments of the computer codes based on the Lanczos procedure for the partial 
eigenanalysis of large symmetric matrices (Ref. 2) motivated similar ongoing developments 
for the unsymmetric case. 

There are two main approaches to this unsymmetric case: the simultaneous iteration 
methods (Refs. 4-7) and the Krylov subspace methods (Refs. 8-11). The success of the 
Lanczos procedures for the symmetric case (Ref. 2) has focused attention on the Krylov 
subspace methods, and it is also increasingly recognized that the simultaneous iteration 
methods are not competitive with the Krylov methods. The Krylov methods comprise two 
sub classes: methods based on orthogonal projection, such as the Arnoldi-Saad or AS 
method (Refs. 8,11), and those based on oblique projection, such as the Lanczos method 
(Ref. 9). Exploring the feasibility of these methods in the Floquet eigenanalysis, even 
within the limited scope of small-order FTMs, should prove useful if only to motivate 
further search for an alternative to the QR method for large rotorcraft systems. The 
present investigation chooses one of the Krylov subspace methods, the AS method, and 
compares it with the QR method for machine-time efficiency and computational reliability. 
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Amoldi-Saad method 


We begin with an n x n FTM, A. The AS method sequentially reduces A to an m x m 
upper Hessenberg matrix H ra ( m << n ), which contains the dominant eigenvalues of A. 
For extracting p dominant eigenpairs, the algorithm generates the Krylov subspace of 
dimension m > p. That is, the algorithm starts with an arbitrary vector vj and generates a 
sequence of vectors vj, V 2 , . .., v m by applying the Gram-Schmidt process to the sequence of 
vectors {v h Avi, A2vi,...,A“*iv,} according to the recurrence formula (Ref. 8): 


hj*t,j v j+ i = Avj - 4 hij vi, 


1=1 


( 1 ) 


The hy (i = l,..,j+l) are chosen so that v j+1 is orthonormal to Vi, i = 1,2,... ,j and 
II v j +i|| = 1- Hence, the set of m vectors {vi, V 2 , •••, v„,} forms an orthornormal basis for m- 
dimensional subspace spanned by { v b Av u A 2 v h ..., A®' 1 vj } . It can be shown (Ref. 8) 
that the n x m matrix B, whose columns comprise m vectors {vi, v 2 , v m }, satisfies 


H m = Bt A B, 


( 2 ) 


where H m is an m x m upper Hessenberg matrix; B will be a tridiagonal matrix if A is 
symmetric. After generating H m , its eigenpairs can be obtained using the QR method. 

Convergence 


In the AS method, the algorithm breaks down at the jth step ( i.e., in the generation of 
the jth Krylov vector) if and only if h j+1 ,j in Eq. (1) becomes zero. But, hj +1 ,j = 0 means 
that the Krylov vectors Ai\ u Ai +1 v, AJ**vj,..., are linearly dependent on the previous j 
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Krylov vectors. The smallest value k of j at which the succeeding Krylov vectors become 
linearly dependent is known as the grade of the vector Vi w.r.t. A (Ref. 12). For such a 
vector of grade k, there are constants a 0 , ai,..,ak such that 

[a 0 I + aiA + a2A 2 + a3A 3 + ,...+akA k ] vi = 0, (3) 

from which it follows that 

[ao + a^ 4- a 2 A 2 + ajA 3 + ,...+akA k ] = 0, (4) 

where A are the eigenvalues of A. Equation (4) is known as the minimal polynomial of vi 
with respect to A (Ref. 12). When Eq. (3) is satisfied, the vectors v b v 2 ,...,vk span the 
dominant subspace of dimension k of A, and the span of those vectors is invariant. Hence, 
the eigenvalues of Hk = B l AB are the exact dominant eigenvalues of A. 

On the other hand, when the number of Krylov vectors chosen is less than the grade of 
the vector v b the accuracy of the eigenvalues computed decreases with their decreasing 
order of magnitude; the most dominant eigenvalue will be the most accurate and the 
smallest eigenvalue will be the least accurate. However, it must be mentioned that issues 
such as roundoff errors and reorthogonalization (Ref. 8) are not addressed in this 
exploratory study. 

Eigenvalue reliability 

We use two eigenvalue reliability parameters: condition number of an eigenvalue and 
residual error of an eigenpair (eigenvalue and corresponding eigenvector). The matrix A 
and its transpose A fc have the same eigenvalues. We assume that A is one such simple 
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eigenvalue (of multiplicity one) and that x and y are the corresponding right and left 
eigenvectors; that is, 


Ax = Ax and A l y = Ay (5) 

Then the condition number of A is given by (Ref. 12) 

cond(A) = | y l x| (6) 

We emphasize that A is usually complex, as are the eigenvectors. It is y l (not the 
hermitian transpose of y) that is used. Further, the vectors are normalized, that is, 
||x|| = || xh x|| = 1 = ||y|| = || y h y||. If A is hermitian, ||yt x|| = 1. The significance of 
cond(A) is demonstrated by the following relation: 

\6X\ < PA|| cond(A), (7) 

where ||&A|| represents the spectral norm of the matrix of perturbations of the FTM due to 
sources such as roundoff and discretization. Further, SX represents the resulting 
perturbation of A due to working with the computed A + £A, instead of A. In other words, 
the computed eigenvalue is A + <5A, which is the exact eigenvalue of A + 6A. Thus, cond(A) 
provides a measure of the sensitivity of an eigenvalue, typified by <5A. If small changes in 
the elements of the FTM can lead to arbitrarily large changes in an eigenvalue, then the 
eigenvalue problem for that eigenvalue is said to be ill-conditioned. That a matrix is well- 
conditioned for eigenvalue computations is no guarantee that it is well-conditioned for 
eigenvector computations. Though the condition number approach has a rigorous basis, a 
similar approach for eigenvectors is too involved for computing the error bounds. 


5 


Therefore, we study the reliability of the computed eigenpair (A,x) by the residual error 
approach, which gives the relative error measure e. That is, 


ll Ax - Ax ll II r II 

e = = ’ W 

II Ax|| || A || 

where ||r|| is the Euclidean norm of the residual error. It appears that cond(A) and e should 
provide adequate information on the reliability of a computed eigenpair. 

Results 

The eigenanalysis comprises computation of a few dominant eigenvalues, including their 
eigenvectors and the corresponding condition numbers and residual errors; see Eqs. (7) and 
(8). These computations from the AS method are compared with those from the QR 
method for FTMs of various order. Table 1 is a representative example. 

The QR results in the last row of Table 1 are taken to be exact; very low residual error 
and condition number are noteworthy. Before comparing the results from these two 
eigenanalysis methods, it is important to realize that the eigenanalysis of large 
unsymmetric matrices represents a computational barrier and several issues concerning the 
minimum dimension of a subspace matrix merit further research. Nevertheless, comparison 
with the QR results shows that the damping and frequency results from the 20 x 20 
subspace matrix are accurate up to six and two significant figures, respectively, with 
comparable condition numbers, and that the accuracy increases with the increasing 
dimension of the subspace matrix. For the AS method, two points are emphasized. First, 
the accuracy is independently corroborated by the respective residual error results. Second, 
the order of magnitude of the eigenvalue condition number is close to the ideal value of 
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one, which demonstrates that the eigenvalue computations are well conditioned. The most 
distinguishing feature is the relative saving in machine time, which is shown in Fig. 1; 53.4 
versus 19.8 seconds for a subspace dimension of 20 x 20. This trend offers great promise for 
large-order FTMs. Because the machine time required for the QR method increases almost 
cubically with the order of the FTM, the relative saving increases rapidly with the 
increasing order. 


Concluding remarks 

For the partial eigenanalyis of the FTM, the AS and the QR methods are compared for 
machine-time efficiency and computational reliability. The numerical experiment shows 
that the Arnoldi-Saad method is more economical than the QR method with comparable 
computational reliability. The partial eigenanalysis of large-order FTMs is important to 
comprehensive stability analysis of rotorcraft, and relatively little is known. Given this 
fact, the results should serve as a useful reference point. 
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Table 1 Machine time and computational reliability results for e 
Arnoldi-Saad method vis-a-vis QR method (FTM dimension - 92x92) 
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Fig. 1 Comparison of machine times required by the 
Arnoldi-Saad and the QR methods (FTM dimension = 92x92) 





























